Self Regulation Ontology DDM Analyses

The tasks included in this report are:

unique(measure_labels$task_group)

 [1] "adaptive_n_back"             "attention_network_task"     
 [3] "choice_reaction_time"        "directed_forgetting"        
 [5] "dot_pattern_expectancy"      "local_global_letter"        
 [7] "motor_selective_stop_signal" "recent_probes"              
 [9] "shape_matching"              "simon"                      
[11] "stim_selective_stop_signal"  "stop_signal"                
[13] "stroop"                      "threebytwo"                 

adaptive_n_back

attention_network_task

choice_reaction_time

directed_forgetting

dot_pattern_expectancy

local_global_letter

motor_selective_stop_signal

recent_probes

shape_matching

simon

stim_selective_stop_signal

stop_signal

stroop

threebytwo

Our analyses span 512 different measures.

NOTE: We found that HDDM parameters for T1 data do not differ in any detectable way depending on whether they are fit using the full T1 sample (n=522) or only for the subset who completed the T2 battery as well (n=150). Therefore the T1 HDDM parameters in this report are based on fits from n=150. Details of these analyses can be found here.

DDM vs raw reliability overall

Plot reliability point estimates comparing DDM measures to raw measures faceting for contrast measures.

fig_name = 'ddmvsraw_point.jpeg'

knitr::include_graphics(paste0(fig_path, fig_name))

Plot averaged bootstrapped reliability estimates per measure comparing DDM measures to raw measures faceting for contrast measures.

fig_name = 'ddmvsraw_boot.jpeg'

knitr::include_graphics(paste0(fig_path, fig_name))

Model testing if the reliability of raw measures differs from that of ddm estimates and if contrast measures differ from non-contrast measures.

Checking if both fixed effects of raw vs ddm and contrast vs non-contrast as well as their interaction is necessary.

Conclusion: Interactive model is best.

mer1 = lmer(icc2.1 ~ ddm_raw + (1|dv), boot_df %>% filter(rt_acc != "other"))
mer1a = lmer(icc2.1 ~ overall_difference + (1|dv), boot_df %>% filter(rt_acc != "other"))
mer2 = lmer(icc2.1 ~ ddm_raw + overall_difference + (1|dv), boot_df %>% filter(rt_acc != "other"))
mer3 = lmer(icc2.1 ~ ddm_raw * overall_difference + (1|dv), boot_df %>% filter(rt_acc != "other"))

fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

anova(mer1, mer2, mer3)

refitting model(s) with ML (instead of REML)

anova(mer1a, mer2, mer3)

refitting model(s) with ML (instead of REML)

rm(mer1, mer2, mer1a)

Raw measures do not significantly differ from ddm parameters in their reliability but non-contrast measures are significantly more reliable compared to contrast and condition measures.

summary(mer3)

Linear mixed model fit by REML ['lmerMod']
Formula: icc2.1 ~ ddm_raw * overall_difference + (1 | dv)
   Data: boot_df %>% filter(rt_acc != "other")

REML criterion at convergence: -546356

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-5.537 -0.630  0.014  0.643  6.045 

Random effects:
 Groups   Name        Variance Std.Dev.
 dv       (Intercept) 0.01841  0.1357  
 Residual             0.00751  0.0866  
Number of obs: 267000, groups:  dv, 267

Fixed effects:
                                      Estimate Std. Error t value
(Intercept)                             0.3681     0.0146   25.16
ddm_rawraw                              0.0837     0.0288    2.91
overall_differencecontrast             -0.2384     0.0218  -10.91
overall_differencecondition            -0.0919     0.0323   -2.85
ddm_rawraw:overall_differencecontrast  -0.0564     0.0397   -1.42

Correlation of Fixed Effects:
                (Intr) ddm_rw ovrll_dffrnccnt ovrll_dffrnccnd
ddm_rawraw      -0.509                                       
ovrll_dffrnccnt -0.670  0.341                                
ovrll_dffrnccnd  0.000 -0.661  0.000                         
ddm_rwrw:v_      0.369 -0.725 -0.550           0.479         
fit warnings:
fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

Best measure for each task

What is the best measure of individual difference for any measure that has both raw and DDM parameters?

Even though overall the ddm parameters are not significantly less reliable the most reliable measure is more frequently a raw measure. There are some examples of an EZ estimate being the best for a task as well. Regardless of raw vs ddm the best measure is always a non-contrast measure.

rel_df %>%
  group_by(task_group) %>%
  filter(icc2.1 == max(icc2.1)) %>%
  select(task_group, everything())

Variance breakdown measure types

fig_name = 'ddmvsraw_varsubs.jpeg'

knitr::include_graphics(paste0(fig_path, fig_name))

fig_name = 'ddmvsraw_varresid.jpeg'

knitr::include_graphics(paste0(fig_path, fig_name))

Model testing if the percentage of between subjects variance of raw measures differs from that of ddm estimates and if contrast measures differ from non-contrast measures.

Checking if both fixed effects of raw vs ddm and contrast vs non-contrast as well as their interaction is necessary.

Conclusion: Model with fixed effects for both is best.

mer1 = lmer(var_subs_pct ~ ddm_raw + (1|dv), boot_df %>% filter(rt_acc != "other"))
mer1a = lmer(var_subs_pct ~ overall_difference + (1|dv), boot_df %>% filter(rt_acc != "other"))
mer2 = lmer(var_subs_pct ~ ddm_raw +  overall_difference + (1|dv), boot_df %>% filter(rt_acc != "other"))
mer3 = lmer(var_subs_pct ~ ddm_raw *  overall_difference + (1|dv), boot_df %>% filter(rt_acc != "other"))

fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

anova(mer1, mer2, mer3)

refitting model(s) with ML (instead of REML)

anova(mer1a, mer2, mer3)

refitting model(s) with ML (instead of REML)

rm(mer1, mer1a, mer3)

Contrast measures have lower between subjects variability (ie are worse individual difference measures). Raw and ddm measures do not differ significantly.

summary(mer2)

Linear mixed model fit by REML ['lmerMod']
Formula: var_subs_pct ~ ddm_raw + overall_difference + (1 | dv)
   Data: boot_df %>% filter(rt_acc != "other")

REML criterion at convergence: 2173929

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.648 -0.725  0.020  0.758  4.603 

Random effects:
 Groups   Name        Variance Std.Dev.
 dv       (Intercept) 142      11.9    
 Residual             200      14.1    
Number of obs: 267000, groups:  dv, 267

Fixed effects:
                            Estimate Std. Error t value
(Intercept)                    42.30       1.20   35.40
ddm_rawraw                      4.15       1.74    2.38
overall_differencecontrast     -6.44       1.60   -4.02
overall_differencecondition     0.53       2.49    0.21

Correlation of Fixed Effects:
                (Intr) ddm_rw ovrll_dffrnccnt
ddm_rawraw      -0.377                       
ovrll_dffrnccnt -0.602 -0.101                
ovrll_dffrnccnd -0.216 -0.519  0.360         

Model testing if the percentage of residual variance of raw measures differs from that of ddm estimates and if contrast measures differ from non-contrast measures.

Checking if both fixed effects of raw vs ddm and contrast vs non-contrast as well as their interaction is necessary.

Conclusion: Interactive model is best

mer1 = lmer(var_resid_pct ~ ddm_raw + (1|dv), boot_df %>% filter(rt_acc != "other"))
mer1a = lmer(var_resid_pct ~ overall_difference + (1|dv), boot_df %>% filter(rt_acc != "other"))
mer2 = lmer(var_resid_pct ~ ddm_raw + overall_difference + (1|dv), boot_df %>% filter(rt_acc != "other"))
mer3 = lmer(var_resid_pct ~ ddm_raw * overall_difference + (1|dv), boot_df %>% filter(rt_acc != "other"))

fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

anova(mer1, mer2, mer3)

refitting model(s) with ML (instead of REML)

anova(mer1a, mer2, mer3)

refitting model(s) with ML (instead of REML)

rm(mer1, mer1a, mer2)

Both contrast and condition measures have higher residual variance. Raw and ddm measures do not differ.

summary(mer3)

Linear mixed model fit by REML ['lmerMod']
Formula: var_resid_pct ~ ddm_raw * overall_difference + (1 | dv)
   Data: boot_df %>% filter(rt_acc != "other")

REML criterion at convergence: 1907553

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-4.042 -0.640 -0.018  0.658  4.753 

Random effects:
 Groups   Name        Variance Std.Dev.
 dv       (Intercept) 49.6     7.04    
 Residual             73.7     8.59    
Number of obs: 267000, groups:  dv, 267

Fixed effects:
                                      Estimate Std. Error t value
(Intercept)                             18.746      0.760   24.67
ddm_rawraw                              -0.973      1.494   -0.65
overall_differencecontrast               9.462      1.134    8.34
overall_differencecondition              3.885      1.676    2.32
ddm_rawraw:overall_differencecontrast    0.501      2.061    0.24

Correlation of Fixed Effects:
                (Intr) ddm_rw ovrll_dffrnccnt ovrll_dffrnccnd
ddm_rawraw      -0.509                                       
ovrll_dffrnccnt -0.670  0.341                                
ovrll_dffrnccnd  0.000 -0.661  0.000                         
ddm_rwrw:v_      0.369 -0.725 -0.550           0.479         
fit warnings:
fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

rm(mer3)

Effects of sample size on reliability can be found here.

Hierarchical estimation consequences

The H in HDDM estimates stands for ‘hierarchical’ denoting the fact that the distribution of the whole sample is incorportated in the priors for the model parameters. This is different than e.g. EZ-diffusion parameters that would be the same for a given subject’s data regardless of the rest of the sample. One can ask whether this H indeed has a measurable impact.

HDDM and EZ parameter estimates might differ due to other difference in the estimation process as well. Therefore to evaluate whether the ‘H’ leads to meaningful changes in parameters we compare the same models fit on the same data using either the whole sample for the hierarchical structure or only the single subject’s data. These latter estimates are referred to as ‘flat’ estimates.

retest_hddm_flat = read.csv(paste0(retest_data_path,'retest_hddm_flat.csv'))

retest_hddm_flat = retest_hddm_flat %>% rename(sub_id = subj_id)

test_hddm_flat = read.csv(paste0(retest_data_path,'/t1_data/t1_hddm_flat.csv'))

test_hddm_flat = test_hddm_flat %>% rename(sub_id = subj_id)

# numeric_cols = get_numeric_cols()

# Check if all the variables are there (no for now)
# sum(names(retest_hddm_flat) %in% numeric_cols) == length(names(retest_hddm_flat))
# names(retest_hddm_flat)[which(names(retest_hddm_flat) %in% numeric_cols == FALSE)]

# sum(names(test_hddm_flat) %in% numeric_cols) == length(names(test_hddm_flat))
# names(test_hddm_flat)[which(names(test_hddm_flat) %in% numeric_cols == FALSE)]

Parameter value

Plot percent change in raw parameters for retest and t1 for each of 3 parameters (using hierarchical as baseline: what percent does the flat parameter change compared to the hierarchical)

For both time points the thresholds and non-decision times are very similary regardless of whether they are estimated heirarchically or without the hierarchy (difference peaking at 0).

This is also true for the majority of the drift rates. There are, however, almost as many drift rates that change completely when estimated without the hierarchy.

#Should not be necessary later
common_cols = names(retest_hddm_flat)[names(retest_hddm_flat) %in% names(retest_data)]
common_cols=common_cols[common_cols %in% names(test_hddm_flat)]

retest_hddm_hier = retest_data %>% select(common_cols) %>% mutate(hddm="hierarchical")
test_hddm_hier = test_data %>% select(common_cols)  %>% mutate(hddm="hierarchical")
retest_hddm_flat = retest_hddm_flat %>% select(common_cols) %>% mutate(hddm="flat")
test_hddm_flat = test_hddm_flat %>% select(common_cols) %>% mutate(hddm="flat")

retest_flat_difference = rbind(retest_hddm_hier, retest_hddm_flat)
retest_flat_difference = retest_flat_difference %>%
  gather(dv, value, -sub_id, -hddm) %>%
  spread(hddm, value) %>%
  mutate(diff_pct = (hierarchical - flat)/hierarchical*100,
         diff_pct = ifelse(diff_pct<(-100), -100, ifelse(diff_pct>100, 100, diff_pct)),
         time = "retest",
         par = ifelse(grepl("drift", dv), "drift", ifelse(grepl("thresh", dv), "thresh", ifelse(grepl("non_decision", dv), "non_decision", NA))))

test_flat_difference = rbind(test_hddm_hier, test_hddm_flat)
test_flat_difference = test_flat_difference %>%
  gather(dv, value, -sub_id, -hddm) %>%
  spread(hddm, value) %>%
  mutate(diff_pct = (hierarchical - flat)/hierarchical*100,
         diff_pct = ifelse(diff_pct<(-100), -100, ifelse(diff_pct>100, 100, diff_pct)),
         time = "test",
         par = ifelse(grepl("drift", dv), "drift", ifelse(grepl("thresh", dv), "thresh", ifelse(grepl("non_decision", dv), "non_decision", NA))))

flat_difference = rbind(test_flat_difference, retest_flat_difference)

flat_difference %>%
  ggplot(aes(diff_pct))+
  geom_histogram()+
  facet_grid(factor(time, levels = c("test", "retest"),labels = c("test", "retest")) ~ par, scales="free")+
  xlab("Percent change")

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: Removed 1022 rows containing non-finite values (stat_bin).

Is the average percentage change different than 0? No.

This model uses the distribution of percentage difference in parameter estimates for drift rates as the baseline. A more appropriate model would test whether all three distributions are different than 0. The conclusion should not change: The intercept suggests that the mean of the drift rate difference distribution is at 0 and this is not different than the distributions for either of the other parameters and either time point.

summary(lmer(diff_pct ~ time*par+(1|dv), flat_difference %>% mutate(pars = factor(par, levels = c("thresh", "drift", "non_decision")))))

Linear mixed model fit by REML ['lmerMod']
Formula: diff_pct ~ time * par + (1 | dv)
   Data: 
flat_difference %>% mutate(pars = factor(par, levels = c("thresh",  
    "drift", "non_decision")))

REML criterion at convergence: 238677

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-3.244 -0.258 -0.016  0.213  3.853 

Random effects:
 Groups   Name        Variance Std.Dev.
 dv       (Intercept)  130     11.4    
 Residual             1443     38.0    
Number of obs: 23578, groups:  dv, 82

Fixed effects:
                         Estimate Std. Error t value
(Intercept)                 -1.19       1.64   -0.73
timetest                     1.64       0.62    2.65
parnon_decision              1.65       3.57    0.46
parthresh                   -1.84       3.39   -0.54
timetest:parnon_decision    -1.38       1.35   -1.02
timetest:parthresh          -1.71       1.28   -1.33

Correlation of Fixed Effects:
            (Intr) timtst prnn_d prthrs tmts:_
timetest    -0.191                            
parnon_dcsn -0.460  0.088                     
parthresh   -0.485  0.093  0.223              
tmtst:prnn_  0.088 -0.459 -0.191 -0.043       
tmtst:prthr  0.092 -0.483 -0.042 -0.192  0.222

Do people change in the same way? Mostly.

Plotting the raw parameter estimate that is estimated hierarchically against the estimate that is estimated without the hierarchy. Red lines are the 45-degree line. The fact that many points cluster around this line and that there are no systematic deviances from it for any parameter at either time point suggests that the hierarchical and flat estimates are mostly the same.

fig_name = 'HDDM_par_flatvshier.jpeg'

knitr::include_graphics(paste0(fig_path, fig_name))

Parameter reliability

Plotting the reliability estimates for flat parameters against the reliability of hierarchical estimates. Red lines are 45-degree lines. While there are some changes in reliability nothing appears very large (or consequential in pushing the reliability to any acceptable level depending on the estimation method) or systematic.

rel_df_flat = make_rel_df(t1_df = test_hddm_flat, t2_df = retest_hddm_flat, metrics = c('icc2.1', 'pearson', 'var_breakdown'))

rel_df_flat = rel_df_flat %>%
  left_join(rel_df[,c("dv", "icc2.1", "rt_acc", "overall_difference")], by = "dv") 

fig_name = 'HDDM_rel_flatvshier.jpeg'

knitr::include_graphics(paste0(fig_path, fig_name))

with(rel_df_flat %>% filter(rt_acc == "drift rate"), t.test(icc2.1.x, icc2.1.y, paired=T))

    Paired t-test

data:  icc2.1.x and icc2.1.y
t = 0.9, df = 51, p-value = 0.4
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.01310  0.03422
sample estimates:
mean of the differences 
                0.01056 

with(rel_df_flat %>% filter(rt_acc == "threshold"), t.test(icc2.1.x, icc2.1.y, paired=T))

    Paired t-test

data:  icc2.1.x and icc2.1.y
t = -1.9, df = 15, p-value = 0.07
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.073416  0.003664
sample estimates:
mean of the differences 
               -0.03488 

with(rel_df_flat %>% filter(rt_acc == "non-decision"), t.test(icc2.1.x, icc2.1.y, paired=T))

    Paired t-test

data:  icc2.1.x and icc2.1.y
t = -0.39, df = 13, p-value = 0.7
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.04242  0.02946
sample estimates:
mean of the differences 
              -0.006479 

rm(flat_difference, rel_df_flat, retest_flat_difference, retest_hddm_flat, retest_hddm_hier, test_flat_difference,  test_hddm_flat, test_hddm_hier)

Parameter fit

Do the fit statistics differ by whether the model was hierarchical or not? No.

t1_hierarchical_fitstats = t1_hierarchical_fitstats %>%
  filter(subj_id %in% retest_hierarchical_fitstats$subj_id)

tmp = rbind(retest_hierarchical_fitstats, retest_flat_fitstats, t1_hierarchical_fitstats, t1_flat_fitstats) %>%
  select(m_kl, subj_id, task_name, sample) %>%
  separate(sample, c("time", "proc"), sep="_", remove=FALSE)

Warning in `[<-.factor`(`*tmp*`, ri, value = c(5L, 7L, 9L, 11L, 12L, 14L, :
invalid factor level, NA generated

Warning in `[<-.factor`(`*tmp*`, ri, value = c(5L, 7L, 9L, 11L, 12L, 14L, :
invalid factor level, NA generated

fig_name = 'HDDM_fitstats_flatvshier.jpeg'

knitr::include_graphics(paste0(fig_path, fig_name))

tmp

summary(lmer(m_kl ~ time*proc+(1|task_name),tmp))

Linear mixed model fit by REML ['lmerMod']
Formula: m_kl ~ time * proc + (1 | task_name)
   Data: tmp

REML criterion at convergence: 30849

Scaled residuals: 
   Min     1Q Median     3Q    Max 
-2.058 -0.429 -0.199  0.163 13.706 

Random effects:
 Groups    Name        Variance Std.Dev.
 task_name (Intercept) 1.11     1.05    
 Residual              2.55     1.60    
Number of obs: 8147, groups:  task_name, 14

Fixed effects:
                        Estimate Std. Error t value
(Intercept)               1.4481     0.2834    5.11
timet1                    0.0487     0.0494    0.99
prochierarchical         -0.1476     0.0498   -2.96
timet1:prochierarchical  -0.0115     0.0710   -0.16

Correlation of Fixed Effects:
            (Intr) timet1 prchrr
timet1      -0.087              
prochrrchcl -0.086  0.495       
tmt1:prchrr  0.061 -0.695 -0.699

Dimensionality reduction

Correlations between variables of same/different type/task

– Do DDM parameters capture similar processes as the raw measures in a given task or do they capture processes that are more similar across tasks? (If the former they would be less useful than if the latter.)

This could be analyzed with factor analysis but there are more variables than observations so as a first pass we’ll explore correlations.

if(!exists('all_data_cor')){
  source('/Users/zeynepenkavi/Dropbox/PoldrackLab/SRO_DDM_Analyses/code/workspace_scripts/var_cor_data.R')
}

fig_name = 'ddm_raw_vars_cor.jpeg'

knitr::include_graphics(paste0(fig_path, fig_name))

(Absolute) correlations between raw and ddm measures within a task are higher than those between ddm measures across tasks. Note that only variables of same type are marked as ddm-ddm or raw-raw (e.g. a correlation between a threshold and another threshold is marked ddm-ddm while a correlation between a threshold and drift rate is not.).

Are the correlations between ddm parameters across tasks higher than correlations between ddm parameters and other variables from a given task? No. Correlations between variables of the same task are consistently higher.

In some sense that is good. We haven’t created all of these different tasks that putatively measure different things for nothing.

summary(lm(log(abs(value)) ~ ddm_ddm*task_task, all_data_cor))

Call:
lm(formula = log(abs(value)) ~ ddm_ddm * task_task, data = all_data_cor)

Residuals:
   Min     1Q Median     3Q    Max 
-9.118 -0.517  0.281  0.794  2.006 

Coefficients:
                                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)                        -2.1047     0.0138 -152.37  < 2e-16 ***
ddm_ddmraw-raw                     -0.2077     0.0196  -10.60  < 2e-16 ***
task_tasksame task                  1.1197     0.0439   25.51  < 2e-16 ***
ddm_ddmraw-raw:task_tasksame task  -0.3610     0.0650   -5.55  2.9e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.11 on 14222 degrees of freedom
  (15784 observations deleted due to missingness)
Multiple R-squared:  0.071, Adjusted R-squared:  0.0708 
F-statistic:  362 on 3 and 14222 DF,  p-value: <2e-16

Factor analysis

Can we recover a 3 factor structure for DDM parameters that captures comparable processes across tasks?

Although the above correlation analyses suggest that variables from the same task are more similar compared to variables across tasks having even one measure per task is not very helpful with numerous tasks (in our case 14) particularly in terms of capturing cognitive mechanisms that are both common and distinct to these tasks.

Therefore we reduce the number of variables by running a factor analysis on a set of variables that ought to capture common processes across the tasks: the DDM parameters.

In another notebook I explore this more extensively and decide to use a 3 factor solution for the EZ variables. To summarize my conclusions:
- A 3 factor solution for the T1 EZ data is not necessarily the model with the lowest BIC but it is not significantly worse from a more complicated model with the lowest BIC.
- HDDM parameters don’t lend themselves as well to this analysis because the majority of the parameters are drift rates since only this parameter was allowed to vary across conditions. Therefore a 3 factor solution on the whole set of HDDM parameters mostly tries to separate out different drift rates from each other.
- Even when using a 3 factor solution the HDDM parameters don’t cluster as separably by parameter types as EZ parameters do.

Here is the depiction of the 3 factor solution for the EZ variables of T1 data.
Bars are colored by the variable type. Bars that are outlined in black load negatively on the factor they are in.

The figure suggests that the three parameter types are largely separable. The first factors captures a large majority of the drift rates. It, however, also contains numerous thresholds that load on negatively. The second factor in the middle contains mostly thresholds and the third factors contains only non-decision times.

fig_name = 'EZ_FA_T1_3.jpeg'

knitr::include_graphics(paste0(fig_path, fig_name))

if(!exists('ez_t1_fa_3')){
  source('/Users/zeynepenkavi/Dropbox/PoldrackLab/SRO_DDM_Analyses/code/workspace_scripts/ez_fa_data.R')
}

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Dropping 0 variables with correlations above 0.85

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Using transformation:log
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
No positively skewed variables found.* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
No negatively skewed variables found.* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Dropping 0 variables with correlations above 0.85

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Using transformation:log
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
No positively skewed variables found.* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
No negatively skewed variables found.* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Dropping 0 variables with correlations above 0.85

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Using transformation:log
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
No positively skewed variables found.* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
No negatively skewed variables found.* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

ez_t1_fa_3 = fa(res_clean_test_data_ez, 3, rotate='oblimin', fm='minres', scores='Anderson')
summary(ez_t1_fa_3)

Factor analysis with Call: fa(r = res_clean_test_data_ez, nfactors = 3, rotate = "oblimin", 
    scores = "Anderson", fm = "minres")

Test of the hypothesis that 3 factors are sufficient.
The degrees of freedom for the model is 21318  and the objective function was  469.9 
The number of observations was  150  with Chi Square =  35790  with prob <  0 

The root mean square of the residuals (RMSA) is  0.09 
The df corrected root mean square of the residuals is  0.09 

Tucker Lewis Index of factoring reliability =  0.858
RMSEA index =  0.124  and the 10 % confidence intervals are  0.066 NA
BIC =  -71027
 With factor correlations of 
      MR1   MR3   MR2
MR1  1.00 -0.20 -0.04
MR3 -0.20  1.00  0.18
MR2 -0.04  0.18  1.00

Reliability of lower dimensions

Are these clusters more reliable than using either the raw or the DDM measures alone?

I had some questions on how to calcualte reliabilities for latent variables. Details of my explorations on these questions can be found here:

I concluded to present two approaches to evaluate this:

1. Predicting T2 using the 3 factor model from the T1 data

fig_name = 'ez_fa_rel_t2pred.jpeg'

knitr::include_graphics(paste0(fig_path, fig_name))

2. Fitting a 3 factor separately for the T2 data and comparing it to the T1 solution

fig_name = 'ez_fa_rel_t2fit.jpeg'

knitr::include_graphics(paste0(fig_path, fig_name))

Here’s how the loadings the two models compare to each other

fig_name = 'ez_fa_t1vst2_loadings.jpeg'

knitr::include_graphics(paste0(fig_path, fig_name))

Prediction

Do raw or DDM measures (or factor scores) predict real world outcomes better?

fig_name = 't1_pred.jpeg'

knitr::include_graphics(paste0(fig_path, fig_name))